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SYSTEM AND METHOD FOR QUANTIFYING TISSUE STRUCTURES AND 

THEIR CHANGE OVER TIME 
Reference to Related Applications 

The present application claims the benefit of U.S. Provisional Application No. 
5 60/306,166, filed July 19, 2001, whose disclosure is hereby incorporated by reference in its 
entirety into the present disclosure. 
Field of the Invention 

The present invention is directed to a system and method for quantifying tissue 
structures and their change over time and is more particularly directed to such a system and 

1 0 method which use biomarkers. 
Description of Related Art 

The measurement of internal organs and structures from CT, MRI, ultrasound, PET, 
and other imaging data sets is an important objective in many fields of medicine. For 
example, in obstetrics, the measurement of the biparietal diameter of the fetal head gives an 

15 objective indicator of fetal growth. Another example is the measurement of the hippocampus 
in patients with epilepsy to determine asymmetry (Ashton E.A., Parker K.J., Berg M.J., and 
Chen C.W. "A Novel Volumetric Feature Extraction Technique with Applications to MR 
Images," IEEE Transactions on Medical Imaging 16:4, 1997). The measurement of the 
thickness of the cartilage of bone is another research area (Stammberger, T., Eckstein, F., 

20 Englmeier, K-H., Reiser, M. "Determination of 3D Cartilage Thickness Data from MR 
Imaging: Computational Method and Reproducibility in the Living," Magnetic Resonance in 
Medicine 41, 1999; and Stammberger, T., Hohe, J., Englmeier, K-H., Reiser, M., Eckstein, F. 
"Elastic Registration of 3D Cartilage Surfaces from MR Image Data for Detecting Local 
Changes in Cartilage Thickness," Magnetic Resonance in Medicine 44, 2000). Those 

25 measurements are quantitative assessments that, when used, are typically based on manual 
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intervention by a trained technician or radiologist. For example, trackball or mouse user 
interfaces are commonly used to derive measurements such as the biparietal diameter. User- 
assisted interfaces are also employed to initiate some semi-automated algorithms (Ashton et 
al). The need for intensive and expert manual intervention is a disadvantage, since the 
5 demarcations can be tedious and prone to a high inter- and intra-observer variability. 
Furthermore, the typical application of manual measurements within 2D slices, or even 
sequential 2D slices within a 3D data-set, is not optimal, since tortuous structures, curved 
structures, and thin structures are not well characterized within a single 2D slice, leading 
again to operator confusion and high variability in results. 

10 The need for accurate and precise measurements of organs, tissues, structures, and 

sub-structures continues to increase. For example, in following the response of a disease to a 
new therapy, the accurate representation of 3D structures is vital in broad areas such as 
neurology, oncology, orthopedics, and urology. Another important need is to track those 
measurements of structures over time, to determine if, for example, a tumor is shrinking or 

15 growing, or if the thin cartilage is further deteriorating. If the structures of interest are 
tortuous, or thin, or curved, or have complicated 3D shapes, then the manual determination of 
the structure from 2D slices is tedious and prone to errors. If those measurements are 
repeated over time on successive scans, then inaccurate trend information can unfortunately 
be obtained. For example, subtle tumor growth along an out-of-plane direction can be lost 

20 within poor accuracy and precision and high variability from manual or semi-manual 
measurements. 

Yet another problem with conventional methods is that they lack sophistication and 
are based on "first order" measurements of diameter, length, or thickness. With some semi- 
manual tracings, the measurement is extended to a two-dimensional area or a three- 
25 dimensional volume (Ashton et al). Those traditional measurements can be insensitive to 



2 



WO 03/009214 PCT/US02/22706 

small but important changes. For example, consider the case of a thin structure such as the 
cartilage. Conventional measurements of volume and thickness will be insensitive to the 
presence or absence of small pits in the cartilage, yet those defects could be an important 
indicator of a disease process. 

5 The prior art is capable of assessing gross abnormalities or gross changes over time. 

However, the conventional measurements are not well suited to assessing and quantifying 
subtle abnormalities, or subtle changes, and are incapable of describing complex topology or 
shape in an accurate manner. Furthermore, manual and semi-manual measurements from raw 
images suffer from a high inter-space and intra-observer variability. Also, manual and semi- 

10 manual measurements tend to produce ragged and irregular boundaries in 3D, when the 
tracings are based on a sequence of 2D images. 
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Summary of the Invention 

It will be readily apparent from the above that a need exists in the art for 
measurements, parameters, and descriptors which are more sophisticated and more 
representative and more sensitive to subtle changes than the simple "first order" 
5 measurements of length, diameter, thickness, area and volume. There is a clear need for 
measurements that are more accurate and precise, with lower variability than conventional 
manual or semi-manual methods. There is furthermore a need for measurements that are 
accurate over time, as repeated measurements are made. There is furthermore a need for 
measurements based on high-resolution data sets, such that small defects, tortuous objects, 

10 thin objects, and curved objects, can be quantified. 

It is therefore a primary object of the invention to provide a more accurate 
quantification of tissue structures. It is another object of the invention to provide a more 
accurate quantification of changes in time of tissue structures. It is a further object of the 
invention to address the needs noted above. 

15 To achieve the above and other objects, the present invention identifies important 

structures or substructures, their normalities and abnormalities, and their specific topological 
and morphological characteristics which are sensitive indicators of joint disease and the state 
of pathology. The abnormality and normality of structures, along with their topological and 
morphological characteristics and radiological and pharmacokinetic parameters, are called 

20 biomarkers, and specific measurements of the biomarkers serve as the quantitative 
assessment of conditions and/or disease. 

In human and animal anatomy texts, there are a great number of named organs, 
structures, and substructures. Furthermore, in disease states modifications to normal 
structures are possible and additional pathological structures or lesions can be present. 

25 Despite the imposing number of defined substructures and pathologies, within the major 
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disease categories and degenerative and other abnormal conditions, there are specific 
parameters that serve as indicators of condition. For example, liver metastases, brain lesions, 
athlerosclerotic plaques, and meniscal tears are some examples of specific indicators of 
different conditions. Those specific indicators are defined as biomarkers of disease. The 
quantification of biomarkers includes the assessment of structural, surface, radiological, and 
pharmacokinetic characteristics that have a non-periodic progression. The accurate and 
sophisticated measurement of biomarkers, and the accurate definition of trends over time, is 
the subject of the present invention. Examples of biomarkers that are measured in the present 
invention are given below. The following lists are intended to be illustrative but not limiting. 
The following biomarkers relate to cancer studies: 

• Tumor surface area 

• Tumor compactness (surface-to-volume ratio) 

• Tumor surface curvature 

• Tumor surface roughness 

• Necrotic core volume 

• necrotic core compactness 

• necrotic core shape 

• Viable periphery volume 

© Volume of tumor vasculature 

• Change in tumor vasculature over time 

• Tumor shape, as defined through spherical harmonic analysis 

• Morphological surface characteristics 

• lesion characteristics 

• tumor characteristics 

• tumor peripheral charact eri sties 
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• tumor core characteristics 

• bone metastases characteristics 

• ascites characteristics 

• pleural fluid characteristics 

5 • vessel structure characteristics 

• neovasculature characteristics 

• polyp characteristics 

• nodule characteristics 

• angiogenisis characteristics 
10 • tumor length 

• tumor width 

• tumor 3D volume 

The following biomarkers are sensitive indicators of osteoarthritis joint disease in 
humans and in animals: 

1 5 • shape of the subchondral bone plate 

• layers of the cartilage and their relative size 

• signal intensity distribution within the cartilage layers 

• contact area between the articulating cartilage surfaces 

• surface topology of the cartilage shape 
20 • intensity of bone marrow edema 

• separation distances between bones 

• meniscus shape 

• meniscus surface area 

• meniscus contact area with cartilage 
25 • cartilage structural characteristics 
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• cartilage surface characteristics 

• meniscus structural characteristics 

• meniscus surface characteristics 

• pannus structural characteristics 
5 • joint fluid characteristics 

• osteophyte characteristics 

• bone characteristics 

• lytic lesion characteristics 

• prosthesis contact characteristics 
10 • prosthesis wear 

• joint spacing characteristics 

• tibia medial cartilage volume 

• Tibia lateral cartilage volume 

• femur cartilage volume 
1 5 • patella cartilage volume 

• tibia medial cartilage curvature 

• tibia lateral cartilage curvature 

• femur cartilage curvature 

• patella cartilage curvature 
20 • cartilage bending energy 

• subchondral bone plate curvature 

• subchondral bone plate bending energy 

• meniscus volume 

• osteophyte volume 

25 • cartilage T2 lesion volumes 
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• bone marrow edema volume and number 

• synovial fluid volume 

• synovial thickening 

• subchondrial bone cyst volume 
5 • kinematic tibial translation 

• kinematic tibial rotation 

• kinematic tibial valcus 

• distance between vertebral bodies 

• degree of subsidence of cage 

1 0 • degree of lordosis by angle measurement 

• degree of off-set between vertebral bodies 

• femoral bone characteristics 

• patella characteristics 

The following new biomarkers are sensitive indicators of neurological disease in 
1 5 humans and in animals: 

• The shape, topology, and morphology of brain lesions 

• The shape, topology, and morphology of brain plaques 

• The shape, topology, and morphology of brain ischemia 

• The shape, topology, and morphology of brain tumors 
20 • The spatial frequency distribution of the sulci and gyri 

• The compactness (a measure of surface to volume ratio) of gray matter and white matter 

• whole brain characteristics 

• gray matter characteristics 
25 • white matter characteristics 

8 
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• cerebral spinal fluid characteristics 

• hippocampus characteristics 

• brain sub-structure characteristics 

• The ratio of cerebral spinal fluid volume to gray mater and white matter volume 
5 • The number and volume of brain lesions 

The following biomarkers are sensitive indicators of disease and toxicity in organs 



• organ volume 

• organ surface 

10 • organ compactness 

• organ shape 

• organ surface roughness 

• fat volume and shape 

Another feature which may be used in the present invention is that of "higher order" 
15 measures. Although the conventional measures of length, diameter, and their extensions to 
area and volume are useful quantities, they are limited in their ability to assess subtle but 
potentially important features of tissue structures or substructures. The example of the 
cartilage was already mentioned, where measures of gross thickness or volume would be 
insensitive to the presence or absence of small defects. Thus, the present invention preferably 
20 uses "higher order" measures of structure and shape to characterize biomarkers. "Higher 
order" measures are defined as any measurements that cannot be extracted directly from the 
data using traditional manual or semi-automated techniques, and that go beyond simple pixel 
counting and that apply directly to 3D and 4D analysis. (Length, area, and volume 
measurements are examples of simple first-order measurements that can be obtained by pixel 
25 counting.) Those higher order measures include, but are not limited to: 
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eigenfunction decompositions 
moments of inertia 

shape analysis, including local curvature 
surface bending energy 
shape signatures 

results of morphological operations such as skeletonization 
fractal analysis 
3D wavelet analysis 

advanced surface and shape analysis such as a 3D orthogonal basis function with scale 
invariant properties 

trajectories of bones, joints, tendons, and moving musculoskeletal structures. 
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Brief Description of the Drawings 

A preferred embodiment of the present invention will be set forth in detail with 
reference to the drawings, in which: 

Fig. 1 shows a flow chart of an overview of the process of the preferred embodiment; 
Fig. 2 shows a flow chart of a segmentation process used in the process of Fig. 1; 
Fig. 3 shows a process of tracking a segmented image in multiple images taken over 

time; 

Fig. 4 shows a block diagram of a system on which the process of Figs. 1-3 can be 
implemented; and 

Fig. 5 shows an image of a biomarker formed in accordance with the preferred 
embodiment. 
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Detailed Description of the Preferred Embodiment 

A preferred embodiment of the present invention will now be set forth with reference 
to the drawings. 

Fig. 1 shows an overview of the process of identifying biomarkers and their trends 
5 over time. In step 102, a three-dimensional image of the organ is taken. In step 104, at least 
one biomarker is identified in the image; the technique for doing so will be explained with 
reference to Fig. 2. In step 106, multiple three-dimensional images of the same region of the 
organ are taken over time. In some cases, step 106 may be completed before step 104; the 
order of the two steps is a matter of convenience. In step 108, the same biomarker or 

10 biomarkers are identified in the images taken over time; the technique for doing so will be 
explained with reference to Fig. 3. The identification of the biomarkers in the multiple image 
allows the development in step 1 10 of a model of the organ in four dimensions, namely, three 
dimensions of space and one of time. From that model, the development of the biomarker or 
biomarkers can be tracked over time in step 112. 

15 The preferred method for extracting the biomarkers is with statistical based reasoning 

as defined in Parker et al (US Patent 6,169,817), whose disclosure is hereby incorporated by 
reference in its entirety into the present disclosure. From raw image data obtained through 
magnetic resonance imaging or the like, an object is reconstructed and visualized in four 
dimensions (both space and time) by first dividing the first image in the sequence of images 

20 into regions through statistical estimation of the mean value and variance of the image data 
and joining of picture elements (voxels) that are sufficiently similar and then extrapolating 
the regions to the remainder of the images by using known motion characteristics of 
components of the image (e.g., spring constants of muscles and tendons) to estimate the rigid 
and deformational motion of each region from image to image. The object and its regions 
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can be rendered and interacted with in a four-dimensional (4D) virtual reality environment, 
the four dimensions being three spatial dimensions and time. 

The segmentation will be explained with reference to Fig. 2. First, at step 201, the 
images in the sequence are taken, as by an MRI. Raw image data are thus obtained. Then, at 
5 step 203, the raw data of the first image in the sequence are input into a computing device. 
Next, for each voxel, the local mean value and region variance of the image data are 
estimated at step 205. The connectivity among the voxels is estimated at step 207 by a 
comparison of the mean values and variances estimated at step 205 to form regions. Once the 
connectivity is estimated, it is determined which regions need to be split, and those regions 

10 are split, at step 209. The accuracy of those regions can be improved still more through the 
segmentation relaxation of step 211. Then, it is determined which regions need to be merged, 
and those regions are merged, at step 213. Again, segmentation relaxation is performed, at 
step 215. Thus, the raw image data are converted into a segmented image, which is the end 
result at step 217. Further details of any of those processes can be found in the above-cited 

1 5 Parker et al patent. 

The creation of a 4D model (in three dimensions of space and one of time) will be 
described with reference to Fig. 3. A motion tracking and estimation algorithm provides the 
information needed to pass the segmented image from one frame to another once the first 
image in the sequence and the completely segmented image derived therefrom as described 

20 above have been input at step 301. The presence of both the rigid and non-rigid components 
should ideally be taken into account in the estimation of the 3D motion. According to the 
present invention, the motion vector of each voxel is estimated after the registration of 
selected feature points in the image. 

To take into consideration the movement of the many structures present in a joint, the 

25 approach of the present invention takes into account the local deformations of soft tissues by 
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using a priori knowledge of the material properties of the different structures found in the 
image segmentation. Such knowledge is input in an appropriate database form at step 303. 
Also, different strategies can be applied to the motion of the rigid structures and to that of the 
soft tissues. Once the selected points have been registered, the motion vector of every voxel 
in the image is computed by interpolating the motion vectors of the selected points. Once the 
motion vector of each voxel has been estimated, the segmentation of the next image in the 
sequence is just the propagation of the segmentation of the former image. That technique is 
repeated until every image in the sequence has been analyzed. 

The definition of time and the order of a sequence can be reversed for the purpose of 
the analysis. For example, in a time series of cancer lesions in the liver, there may be more 
lesions in the final scan than were present in the initial scan. Thus, the 4D model can be run 
in the reverse direction to make sure all lesions are accounted for. Similarly, a long time 
series can be run from a mid-point, with analysis proceeding both forward and backward 
from the mid-point. 

Finite-element models (FEM) are known for the analysis of images and for time- 
evolution analysis. The present invention follows a similar approach and recovers the point 
correspondence by minimizing the total energy of a mesh of masses and springs that models 
the physical properties of the anatomy. In the present invention, the mesh is not constrained 
by a single structure in the image, but instead is free to model the whole volumetric image, in 
which topological properties are supplied by the first segmented image and the physical 
properties are supplied by the a priori properties and the first segmented image. The motion 
estimation approach is an FEM-based point correspondence recovery algorithm between two 
consecutive images in the sequence. Each node in the mesh is an automatically selected 
feature point of the image sought to be tracked, and the spring stiffness is computed from the 
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first segmented image and a priori knowledge of the human anatomy and typical 
biomechanical properties for muscle, bone and the like. 

Many deformable models assume that a vector force field that drives spring-attached 
point masses can be extracted from the image. Most such models use that approach to build 
5 semi-automatic feature extraction algorithms. The present invention employs a similar 
approach and assumes that the image sampled at t = n is a set of three dynamic scalar fields: 
<D(x,0=te n W > |Vg n (x)| ) V 2 gn (x)}, 

namely, the gray-scale image value, the magnitude of the gradient of the image value, and the 
Laplacian of the image value. Accordingly, a change in <D(x, t) causes a quadratic change in 

10 the scalar field energy £/<d(x) °c (AO(x)) 2 . Furthermore, the structures underlying the image 
are assumed to be modeled as a mesh of spring-attached point masses in a state of 
equilibrium with those scalar fields. Although equilibrium assumes that there is an external 
force field, the shape of the force field is not important. The distribution of the point masses 
is assumed to change in time, and the total energy change in a time period A/ after time / = n 

15 is given by 
At/ n (Ax) = 

I [(*(*.(*)-*„,(* + Ax))) 2 +(A|V^(x)|~|V g/l+1 (x + Ax)|)) 2 + 
{7iS 2 g n W + V 2 g n+1 (x + Ax))) 2 + ±r,AX T KAX] 

where a, P, and y are weights for the contribution of every individual field change, rj weighs 
the gain in the strain energy, K is the FEM stiffness matrix, and AX is the FEM node 
displacement matrix. Analysis of that equation shows that any change in the image fields or 
20 in the mesh point distribution increases the system total energy. Therefore, the point 
correspondence from g n to g„+i is given by the mesh configuration whose total energy 
variation is a minimum. Accordingly, the point correspondence is given by 

15 
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X = X + AX 
where 

AX = min^ AU n (AX). 

In that notation, min,, q is the value of p that minimizes q. 
5 While the equations set forth above could conceivably be used to estimate the motion 

(point correspondence) of every voxel in the image, the number of voxels, which is typically 
over one million, and the complex nature of the equations make global minimization difficult. 
To simplify the problem, a coarse FEM mesh is constructed with selected points from the 
image at step 305. The energy minimization gives the point correspondence of the selected 
10 points. 

The selection of such points is not trivial. First, for practical purposes, the number of 
points has to be very small, typically = 10 4 ; care must be taken that the selected points 
describe the whole image motion. Second, region boundaries are important features because 
boundary tracking is enough for accurate region motion description. Third, at region 

15 boundaries, the magnitude of the gradient is high, and the Laplacian is at a zero crossing 
point, making region boundaries easy features to track. Accordingly, segmented boundary 
points are selected in the construction of the FEM. 

Although the boundary points represent a small subset of the image points, there are 
still too many boundary points for practical purposes. In order to reduce the number of 

20 points, constrained random sampling of the boundary points is used for the point extraction 
step. The constraint consists of avoiding the selection of a point too close to the points 
already selected. That constraint allows a more uniform selection of the points across the 
boundaries. Finally, to reduce the motion estimation error at points internal to each region, a 
few more points of the image are randomly selected using the same distance constraint. 

25 Experimental results show that between 5,000 and 10,000 points are enough to estimate and 

16 
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describe the motion of a typical volumetric image of 256x256x34 voxels. Of the selected 
points, 75% are arbitrarily chosen as boundary points, while the remaining 25% are interior 
points. Of course, other percentages can be used where appropriate. 

Once a set of points to track is selected, the next step is to construct an FEM mesh for 
5 those points at step 307. The mesh constrains the kind of motion allowed by coding the 
material properties and the interaction properties for each region. The first step is to find, for 
every nodal point, the neighboring nodal point. Those skilled in the art will appreciate that 
the operation of finding the neighboring nodal point corresponds to building the Voronoi 
diagram of the mesh. Its dual, the Delaunay triangulation, represents the best possible 

10 tetrahedral finite element for a given nodal configuration. The Voronoi diagram is 
constructed by a dilation approach. Under that approach, each nodal point in the discrete 
volume is dilated. Such dilation achieves two purposes. First, it is tested when one dilated 
point contacts another, so that neighboring points can be identified. Second, every voxel can 
be associated with a point of the mesh. 

15 Once every point x, has been associated with a neighboring point Xy, the two points are 

considered to be attached by a spring having spring constant kjj 9 where / and m identify the 

materials. The spring constant is defined by the material interaction properties of the 
connected points; those material interaction properties are predefined by the user in 
accordance with known properties of the materials. If the connected points belong to the 
20 same region, the spring constant reduces to kjj and is derived from the elastic properties of 

the material in the region. If the connected points belong to different regions, the spring 
constant is derived from the average interaction force between the materials at the boundary. 
If the object being imaged is a human shoulder, the spring constant can be derived from a 
table such as the following: 

i i i i i l 
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Humeral head 


Muscle 


Tendon 


Cartilage 


Humeral head 


10 4 


0.15 


0.7 


0.01 


Muscle 


0.15 


0.1 


0.7 


0.6 


Tendon 


0.7 


0.7 


10 


0.01 


Cartilage 


0.01 


0.6 


0.01 


10 2 



In theory, the interaction must be defined between any two adjacent regions. In 
practice, however, it is an acceptable approximation to define the interaction only between 
major anatomical components in the image and to leave the rest as arbitrary constants. In 
5 such an approximation, the error introduced is not significant compared with other errors 
introduced in the assumptions set forth above. 

Spring constants can be assigned automatically, as the approximate size and image 
intensity for the bones are usually known a priori. Segmented image regions matching the a 
priori expectations are assigned to the relatively rigid elastic constants for bone. Soft tissues 
10 and growing or shrinking lesions are assigned relatively soft elastic constants. 

Once the mesh has been set up, the next image in the sequence is input at step 309, 
and the energy between the two successive images in the sequence is minimized at step 311. 
The problem of minimizing the energy U can be split into two separate problems: 
minimizing the energy associated with rigid motion and minimizing that associated with 
1 5 deformable motion. While both energies use the same energy function, they rely on different 
strategies. 

The rigid motion estimation relies on the fact that the contribution of rigid motion to 
the mesh deformation energy (AX r KAX)/2 is very close to zero. The segmentation and the a 
priori knowledge of the anatomy indicate which points belong to a rigid body. If such points 

18 
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are selected for every individual rigid region, the rigid motion energy minimization is 
accomplished by finding, for each rigid region the rigid motion rotation R/ and the 
translation T, that minimize that region's own energy: 
= *!< U rigid = £(A* = min^ t/„(A* f )) 

VUrigid 

5 where AX,- = R,X/ + T,Xi and Ax, is the optimum displacement matrix for the points that 

belong to the rigid region R im That minimization problem has only six degrees of freedom for 
each rigid region: three in the rotation matrix and three in the translation matrix. Therefore, 
the twelve components (nine rotational and three translational) can be found via a six- 
dimensional steepest-descent technique if the difference between any two images in the 

10 sequence is small enough. 

Once the rigid motion parameters have been found, the deformational motion is 
estimated through minimization of the total system energy U. That minimization cannot be 
simplified as much as the minimization of the rigid energy, and without further 
considerations, the number of degrees of freedom in a 3D deformable object is three times the 

15 number of node points in the entire mesh. The nature of the problem allows the use of a 
simple gradient descent technique for each node in the mesh. From the potential and kinetic 
energies, the Lagrangian (or kinetic potential, defined in physics as the kinetic energy minus 
the potential energy) of the system can be used to derive the Euler-Lagrange equations for 
every node of the system where the driving local force is just the gradient of the energy field. 

20 For every node in the mesh, the local energy is given by 

(a(g n ( Xi + Ax)-g fl+1 (x,))) 2 +{/3<\VgM + Ax)|-| V ^-i(^)|)> 2 + 

where G m represents a neighborhood in the Voronoi diagram. 
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Thus, for every node, there is a problem in three degrees of freedom whose 
minimization is performed using a simple gradient descent technique that iteratively reduces 
the local node energy. The local node gradient descent equation is 
x t (n + 1) = *,(«) - v&U Mn)in) (Ax) 

5 where the gradient of the mesh energy is analytically computable, the gradient of the field 
energy is numerically estimated from the image at two different resolutions, x(«+l) is the 
next node position, and v is a weighting factor for the gradient contribution. 

At every step in the minimization, the process for each node takes into account the 
neighboring nodes' former displacement. The process is repeated until the total energy 
10 reaches a local minimum, which for small deformations is close to or equal to the global 
minimum. The displacement vector thus found represents the estimated motion at the node 
points. 

Once the minimization process just described yields the sampled displacement field 
AX, that displacement field is used to estimate the dense motion field needed to track the 
15 segmentation from one image in the sequence to the next (step 313). The dense motion is 
estimated by weighting the contribution of every neighbor mode in the mesh. A constant 
velocity model is assumed, and the estimated velocity of a voxel x at a time / is v(x, t) = 
Ax(/)/Af . The dense motion field is estimated by 
k l ' m Ax : 



/ \ c(x) 



& VA*, € <7 W (* ( ) J*"*, | 

20 where 



kf' m is the spring constant or stiffiiess between the materials / and m associated with the voxels 
x and x y , A/ is the time interval between successive images in the sequence, |x - Xy| is the 
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simple Euclidean distance between the voxels, and the interpolation is performed using the 
neighbor nodes of the closest node to the voxel x. That interpolation weights the contribution 
of every neighbor node by its material property kj" ; thus, the estimated voxel motion is 

similar for every homogeneous region, even at the boundary of that region. 
5 Then, at step 315, the next image in the sequence is filled with the segmentation data. 

That means that the regions determined in one image are carried over into the next image. To 
do so, the velocity is estimated for every voxel in that next image. That is accomplished by a 
reverse mapping of the estimated motion, which is given by 

v(*,/ + A0 = -ir £>(* y .,0 

10 where H is the number of points that fall into the same voxel space S(x) in the next image. 
That mapping does not fill all the space at time t+At, but a simple interpolation between 
mapped neighbor voxels can be used to fill out that space. Once the velocity is estimated for 
every voxel in the next image, the segmentation of that image is simply 
L(x, t + AO = L(x - v(x, t + AO At, t) 

15 where L(x 9 t) and Z,(x,/+A/) are the segmentation labels at the voxel x for the times / and /+Af. 

At step 317, the segmentation thus developed is adjusted through relaxation labeling, 
such as that done at steps 211 and 215, and fine adjustments are made to the mesh nodes in 
the image. Then, the next image is input at step 309, unless it is determined at step 319 that 
the last image in the sequence has been segmented, in which case the operation ends at step 

20 321. 

The operations described above can be implemented in a system such as that shown in 
the block diagram of Fig. 4. System 400 includes an input device 402 for input of the image 
data, the database of material properties, and the like. The information input through the 
input device 402 is received in the workstation 404, which has a storage device 406 such as a 
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hard drive, a processing unit 408 for performing the processing disclosed above to provide 
the 4D data, and a graphics rendering engine 410 for preparing the 4D data for viewing, e.g., 
by surface rendering. An output device 412 can include a monitor for viewing the images 
rendered by the rendering engine 410, a further storage device such as a video recorder for 
5 recording the images, or both. Illustrative examples of the workstation 304 and the graphics 
rendering engine 410 are a Silicon Graphics Indigo workstation and an Irix Explorer 3D 
graphics engine. 

Shape and topology of the identified biomarkers can be quantified by any suitable 
techniques known in analytical geometry. The preferred method for quantifying shape and 
10 topology is with the morphological and topological formulas as defined by the following 
references: 

Shape Analysis and Classification, L. Costa and R. Cesar, Jr., CRC Press, 2001. 
Curvature Analysis: Peet, F.G., Sahota, T.S. "Surface Curvature as a Measure of 
Image Texture" IEEE Transactions on Pattern Analysis and Machine Intelligence 1985 Vol 
15 PAMI-7 G:734-738. 

Struik, D.J., Lectures on Classical Differential Geometry, 2nd ed., Dover, 1988. 
Shape and Topological Descriptors: Duda, R.O, Hart, P.R, Pattern Classification and 
Scene Analysis, Wiley & Sons, 1973. 

Jain, A.K, "Fundamentals of Digital Image Processing," Prentice Hall, 1989. 
20 Spherical Harmonics: Matheny, A., Goldgof, D. "The Use of Three and Four 

Dimensional Surface Harmonics for Nonrigid Shape Recovery and Representation," IEEE 
Transactions on Pattern Analysis and Machine Intelligence 1995, 17: 967-981.;. Chen, C.W, 
Huang, T.S., Arrot, M. "Modeling, Analysis, and Visualization of Left Ventricle Shape and 
Motion by Hierarchical Decomposition," IEEE Transactions on Pattern Analysis and 
25 Machine Intelligence 1994, 342-356. 
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Those morphological and topological measurements have not in the past been applied 
to biomarkers which have a progressive, non-periodic change over time. 

As one example of the quantitative measurement of new biomarkers, the knee of an 
adult human was scanned with a 1 .STesla MRI system, with an in-plane resolution of 0.3 mm 
5 and a slice thickness of 2.0 mm. The cartilage of the femur, tibia, and fibia were segmented 
using the statistical reasoning techniques of Parker et al (cited above). Characterization of 
the cartilage structures was obtained by applying morphological and topological 
measurements. One such measurement is the estimation of local surface curvature. 
Techniques for the determination of local surface curvature are well known in analytic 

10 geometry. For example, if S(x,y,z) is the surface of a structure with an outward normal <n> 
the mean curvature, a local quantity can be determined from the roots of a quadratic equation 
found in Struik (cited above), p. 83. The measurements provide a quantitative, reproducible, 
and very sensitive characterization of the cartilage, in a way which is not practical using 
conventional manual tracings of 2D image slices. 

15 Figure 5 provides a gray scale graph of the quantitative higher order measure of 

surface curvature, as a function of location within the surface of the cartilage. The view is 
from the upper femur, looking down towards the knee to the inner surface of the cartilage. 
Shades of dark-to-light indicate quantitative measurements of local curvature, a higher order 
measurement. 

20 Those data are then analyzed over time as the individual is scanned at later intervals. 

There are two types of presentations of the time trends that are preferred. In one class, the 
repeated higher order measurements are as shown as in Fig. 5, with successive measurements 
overlaid in rapid sequence so as to form a movie. In the complementary representation, a 
trend plot is drawn giving the higher order measures as a function of time. For example, the 
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mean and standard deviation (or range) of the local curvature can be plotted for a specific 
cartilage local area, as a function of time. 

The accuracy of those measurements and their sensitivity to subtle changes in small 
substructures are highly dependent on the resolution of the imaging system. Unfortunately, 
5 most CT, MRI, and ultrasound systems have poor resolution in the out-of-plane, or "z" axis. 
While the in-plane resolution of those systems can commonly resolve objects that are just 
under one millimeter in separation, the out-of-plane (slice thickness) is commonly set at 
1.5mm or even greater. For assessing subtle changes and small defects using higher order 
structural measurements, it is desirable to have better than one millimeter resolution in all 

10 three orthogonal axes. That can be accomplished by fusion of a high resolution scan in the 
orthogonal, or out-of-plane direction, to create a high resolution voxel data set (Pefia, J.-T., 
Totterman, S.M.S., Parker, KJ. "MRI Isotropic Resolution Reconstruction from Two 
Orthogonal Scans," SPIE Medical Imaging, 2001, hereby incorporated by reference in its 
entirety into the present disclosure). In addition to the assessment of subtle defects in 

15 structures, that high-resolution voxel data set enables more accurate measurement of 
structures that are thin, curved, or tortuous. 

In following the response of a person or animal to therapy, or to monitor the 
progression of disease, it is desirable to accurately and precisely monitor the trends in 
biomarkers over time. That is difficult to do in conventional practice since repeated scans 

20 must be reviewed independently and the biomarkers of interest must be traced or measured 
manually or semi-manually with each time interval representing a new and tedious process 
for repeating the measurements. It is highly advantageous to take a 4D approach, such as was 
defined in the above-cited patent to Parker et al, where a biomarker is identified with 
statistical reasoning, and the biomarker is tracked from scan to scan over time. That is, the 

25 initial segmentation of the biomarker of interest is passed on to the data sets from scans taken 
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at later intervals. A search is done to track the biomarker boundaries from one scan to the 
next. The accuracy and precision and reproducibility of that approach is superior to that of 
performing manual or semi-manual measurements on images with no automatic tracking or 
passing of boundary information from one scan interval to subsequent scans. 
5 The quantitative assessment of the new biomarkers listed above provides an objective 

measurement of the state of the joints, particularly in the progression of joint disease. It is 
also very useful to obtain accurate measurements of those biomarkers over time, particularly 
to judge the degree of response to a new therapy, or to assess the trends with increasing age. 
Manual and semi-manual tracings of conventional biomarkers (such as the simple thickness 

10 or volume of the cartilage) have a high inherent variability, so as successive scans are traced 
the variability can hide subtle trends. That means that only gross changes, sometimes over 
very long time periods, can be verified in conventional methods. The inventors have 
discovered that by extracting the biomarker using statistical tests, and by treating the 
biomarker over time as a 4D object, with an automatic passing of boundaries from one time 

15 interval to the next, provides a highly accurate and reproducible segmentation from which 
trends over time can be detected. Thus, the combination of selected biomarkers that 
themselves capture subtle pathologies, with a 4D approach to increase accuracy and 
reliability over time, creates sensitivity that has not been previously obtainable. 

While a preferred embodiment of the invention has been set forth above, those skilled 

20 in the art who have reviewed the present disclosure will readily appreciate that other 
embodiments can be realized within the scope of the present invention. For example, any 
suitable imaging technology can be used. Therefore, the present invention should be 
construed as limited only by the appended claims. 
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We claim: 

1 . A method for assessing a region of interest of a patient, the method comprising: 

(a) taking at least one three-dimensional image of the region of interest; 

(b) identifying at least one biomarker in the at least one three-dimensional image; and 

(c) storing the at least one three-dimensional image and an identification of the at least 
one biomarker in a storage medium. 

2. The method of claim 1, wherein step (b) comprises statistical segmentation of the at 
least one three-dimensional image to identify the at least one biomarker. 

3. The method of claim 1, wherein the at least one three-dimensional image comprises 
a plurality of three-dimensional images of the region of interest taken over time. 

4. The method of claim 3, wherein step (b) comprises statistical segmentation of a 
three-dimensional image selected from the plurality of three-dimensional images to identify 
the at least one biomarker. 

5. The method of claim 4, wherein step (b) further comprises motion tracking and 
estimation to identify the at least one biomarker in the plurality of three-dimensional images 
in accordance with the at least one biomarker identified in the selected three-dimensional 
image. 

6. The method of claim 5, wherein the plurality of three-dimensional images and the 
at least one biomarker identified in the plurality of three-dimensional images are used to form 
a model of the region of interest and the at least one biomarker in three dimensions of space 
and one dimension of time. 

7. The method of claim 6, wherein the biomarker is tracked over time in the model. 

8. The method of claim 1, wherein a resolution in all three dimensions of the at least 
one three-dimensional image is finer than 1 mm. 
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9. The method of claim 1, further comprising deriving a quantitative measure of the at 
least one biomarker. 

10. The method of claim 9, wherein the quantitative measure comprises a 
morphological and topological measure. 

11. The method of claim 10, wherein the morphological and topological measurement 
comprises an estimate of local surface curvature. 

12. The method of claim 1, wherein the at least one biomarker is selected from the 
group consisting of: 

tumor surface area; 

tumor compactness; 

tumor surface curvature; 

tumor surface roughness; 

necrotic core volume; 

necrotic core compactness; 

necrotic core shape; 

viable periphery volume; 

volume of tumor vasculature; 

change in tumor vasculature over time; 

tumor shape; 

morphological surface characteristics; 

lesion characteristics; 

tumor characteristics; 

tumor peripheral characteristics; 

tumor core characteristics; 

bone metastases characteristics; 
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ascites characteristics; 

pleural fluid characteristics; 

vessel structure characteristics; 

neovasculature characteristics; 

polyp characteristics; 

nodule characteristics; 

angiogenisis characteristics; 

tumor length; 

tumor width; 

tumor 3d volume; 

shape of a subchondral bone plate; 

layers of cartilage and relative size of said layers; 

signal intensity distribution within cartilage layers; 

contact area between articulating cartilage surfaces; 

surface topology of cartilage shape; 

intensity of bone marrow edema; 

separation distances between bones; 

meniscus shape; 

meniscus surface area; 

meniscus contact area with cartilage; 

cartilage structural characteristics; 

cartilage surface characteristics; 

meniscus structural characteristics; 

meniscus surface characteristics; 

pannus structural characteristics; 
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joint fluid characteristics; 

osteophyte characteristics; 

bone characteristics; 

lytic lesion characteristics; 

prosthesis contact characteristics; 

prosthesis wear; 

joint spacing characteristics; 

tibia medial cartilage volume; 

tibia lateral cartilage volume; 

femur cartilage volume; 

patella cartilage volume; 

tibia medial cartilage curvature; 

tibia lateral cartilage curvature; 

femur cartilage curvature; 

patella cartilage curvature; 

cartilage bending energy; 

subchondral bone plate curvature; 

subchondral bone plate bending energy; 

meniscus volume; 

osteophyte volume; 

cartilage t2 lesion volumes; 

bone marrow edema volume and number; 

synovial fluid volume; 

synovial thickening; 

subchondrial bone cyst volume; 
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kinematic tibial translation; 

kinematic tibial rotation; 

kinematic tibial valcus; 

distance between vertebral bodies; 

degree of subsidence of cage; 

degree of lordosis by angle measurement; 

degree of off-set between vertebral bodies; 

femoral bone characteristics; 

patella characteristics; 

a shape, topology, and morphology of brain lesions; 
a shape, topology, and morphology of brain plaques; 
a shape, topology, and morphology of brain ischemia; 
a shape, topology, and morphology of brain tumors 
a spatial frequency distribution of the sulci and gyri; 
compactness of gray matter and white matter; 
whole brain characteristics; 
gray matter characteristics; 
white matter characteristics; 
cerebral spinal fluid characteristics; 
hippocampus characteristics; 
brain sub-structure characteristics; 

a ratio of cerebral spinal fluid volume to gray mater and white matter volume; 
the number and volume of brain lesions; 
organ volume; 
organ surface; 
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organ compactness; 
organ shape; 

organ surface roughness; and 
fat volume and shape. 

5 13. The method of claim 1, wherein step (b) comprises taking a higher-order measure 

of the at least one biomarker. 

14. The method of claim 13, wherein the higher-order measure is selected from the 
group consisting of: 

eigenfunction decompositions; 
1 0 moments of inertia; 

shape analysis; 
surface bending energy; 
shape signatures; 

results of morphological operations; 
15 fractal analysis; 

3D wavelet analysis; 

advanced surface and shape analysis; and 

trajectories of bones, joints, tendons, and moving musculoskeletal structures. 

15. The method of claim 13, wherein the higher order measure is obtained as a 
20 function of time from a plurality of three-dimensional images. 

16. The method of claim 1, wherein step (a) is performed through magnetic resonance 
imaging. 

17. A system for assessing a region of interest of a patient, the system comprising: 

(a) an input device for receiving at least one three-dimensional image of the region of 
25 interest; 
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(b) a processor, in communication with the input device, for receiving the at least one 
three-dimensional image of the region of interest from the input device and for identifying at 
least one biomarker in the at least one three-dimensional image; 

(c) storage, in communication with the processor, for storing the at least one three- 
5 dimensional image and an identification of the at least one biomarker; and 

(d) an output device for displaying the at least one three-dimensional image and the 
identification of the at least one biomarker. 

18. The system of claim 17, wherein the processor identifies the at least one 

biomarker through statistical segmentation of the at least one three-dimensional image. 
10 19. The system of claim 17, wherein the at least one three-dimensional image 

comprises a plurality of three-dimensional images of the region of interest taken over time. 

20. The system of claim 19, wherein the processor identifies the at least one 

biomarker through statistical segmentation of a three-dimensional image selected from the 

plurality of three-dimensional images. 
15 21. The system of claim 20, wherein the processor uses motion tracking and 

estimation to identify the at least one biomarker in the plurality of three-dimensional images 

in accordance with the at least one biomarker identified in the selected three-dimensional 

image. 

22. The system of claim 21, wherein the plurality of three-dimensional images and the 
20 at least one biomarker identified in the plurality of three-dimensional images are used to form 

a model of the region of interest and the at least one biomarker in three dimensions of space 
and one dimension of time. 

23. The system of claim 17, wherein a resolution in all three dimensions of the at least 
one three-dimensional image is finer than 1 mm. 
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24. The system of claim 17, wherein the processor derives a quantitative measure of 
the at least one biomarker. 

25. The system of claim 24, wherein the quantitative measure comprises a 
morphological and topological measure. 

26. The system of claim 25, wherein the morphological and topological measurement 
comprises an estimate of local surface curvature. 

27. The system of claim 17, wherein the at least one biomarker is selected from the 
group consisting of: 

tumor surface area; 

tumor compactness; 

tumor surface curvature; 

tumor surface roughness; 

necrotic core volume; 

necrotic core compactness; 

necrotic core shape; 

viable periphery volume; 

volume of tumor vasculature; 

change in tumor vasculature over time; 

tumor shape; 

morphological surface characteristics; 

lesion characteristics; 

tumor characteristics; 

tumor peripheral characteristics; 

tumor core characteristics; 

bone metastases characteristics; 
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ascites characteristics; 

pleural fluid characteristics; 

vessel structure characteristics; 

neo vasculature characteristics; 

polyp characteristics; 

nodule characteristics; 

angiogenisis characteristics; 

tumor length; 

tumor width; 

tumor 3d volume; 

shape of a subchondral bone plate; 

layers of cartilage and relative size of said layers; 

signal intensity distribution within cartilage layers; 

contact area between articulating cartilage surfaces; 

surface topology of cartilage shape; 

intensity of bone marrow edema; 

separation distances between bones; 

meniscus shape; 

meniscus surface area; 

meniscus contact area with cartilage; 

cartilage structural characteristics; 

cartilage surface characteristics; 

meniscus structural characteristics; 

meniscus surface characteristics; 

pannus structural characteristics; 
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joint fluid characteristics; 

osteophyte characteristics; 

bone characteristics; 

lytic lesion characteristics; 

prosthesis contact characteristics; 

prosthesis wear; 

joint spacing characteristics; 

tibia medial cartilage volume; 

tibia lateral cartilage volume; 

femur cartilage volume; 

patella cartilage volume; 

tibia medial cartilage curvature; 

tibia lateral cartilage curvature; 

femur cartilage curvature; 

patella cartilage curvature; 

cartilage bending energy; 

subchondral bone plate curvature; 

subchondral bone plate bending energy; 

meniscus volume; 

osteophyte volume; 

cartilage t2 lesion volumes; 

bone marrow edema volume and number; 

synovial fluid volume; 

synovial thickening; 

subchondrial bone cyst volume; 
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kinematic tibial translation; 

kinematic tibial rotation; 

kinematic tibial valcus; 

distance between vertebral bodies; 

degree of subsidence of cage; 

degree of lordosis by angle measurement; 

degree of off-set between vertebral bodies; 

femoral bone characteristics; 

patella characteristics; 

a shape, topology, and morphology of brain lesions; 
a shape, topology, and morphology of brain plaques; 
a shape, topology, and morphology of brain ischemia; 
a shape, topology, and morphology of brain tumors 
a spatial frequency distribution of the sulci and gyri; 
compactness of gray matter and white matter; 
whole brain characteristics; 
gray matter characteristics; 
white matter characteristics; 
cerebral spinal fluid characteristics; 
hippocampus characteristics; 
brain sub-structure characteristics; 

a ratio of cerebral spinal fluid volume to gray mater and white matter volume; 
the number and volume of brain lesions; 
organ volume; 
organ surface; 
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organ compactness; 
organ shape; 

organ surface roughness; and 
fat volume and shape. 

28. The system of claim 17, wherein the processor takes a higher-order measure of the 
at least one biomarker. 

29. The system of claim 28, wherein the higher-order measure is selected from the 
group consisting of: 

eigenfiinction decompositions; 
moments of inertia; 
shape analysis; 
surface bending energy; 
shape signatures; 

results of morphological operations; 

fractal analysis; 

3D wavelet analysis; 

advanced surface and shape analysis; and 

trajectories of bones, joints, tendons, and moving musculoskeletal structures. 

30. The system of claim 28, wherein the higher order measure is obtained as a 
function of time from a plurality of three-dimensional images. 
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□ BLACK BORDERS 

□ IMAGE CUT OFF AT TOP, BOTTOM OR SIDES 



hi FADED TEXT OR DRAWING 

□ BLURRED OR ILLEGIBLE TEXT OR DRAWING 

□ SKEWED/SLANTED IMAGES 

□ COLOR OR BLACK AND WHITE PHOTOGRAPHS 

□ GRAY SCALE DOCUMENTS . 

□ LINES OR MARKS ON ORIGINAL DOCUMENT 

□ REFERENCE(S) OR EXHIBIT(S) SUBMITTED ARE POOR QUALITY 

□ OTHER: 

IMAGES ARE BEST AVAILABLE COPY. 
As rescanning these documents will not correct the image 
problems checked, please do not report these problems to 
the IFW Image Problem Mailbox. 
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